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Abstract A class of high-order numerical algorithms for Riesz derivatives are 
established through constructing new generating functions. Such new high- 
order formulas can be regarded as the modification of the classical (or shifted) 
Lubich’s difference ones, which greatly improve the convergence orders and sta¬ 
bility for time-dependent problems with Riesz derivatives. In rapid sequence, 
we apply the 2nd-order formula to one-dimension Riesz spatial fractional par¬ 
tial differential equations to establish an unconditionally stable finite difference 
scheme with convergent order 0{t 2 + h 2 ), where r and h are the temporal and 
spatial stepsizes, respectively. Finally, some numerical experiments are per¬ 
formed to confirm the theoretical results and testify the effectiveness of the 
derived numerical algorithms. 
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1 Introduction 

In recent years, increasing attentions have been attracted on fractional calcu¬ 
lus due to its widespread applications in science and engineering [1511151 . In 


The work was partially supported by the National Natural Science Foundation of China un¬ 
der Grant Nos. 11372170 and 11561060, the Scientific Research Program for Young Teachers 
of Tianshui Normal University under Grant No. TSA1405, and Tianshui Normal University 
Key Construction Subject Project (Big data processing in dynamic image). 


Hengfei Ding 

School of Mathematics and Statistics, Tianshui Normal University, Tianshui 741001, China 
E-mail: dinghf05@163.com 

Changpin Li 

Department of Mathematics, Shanghai University, Shanghai 200444, China 
E-mail: lcp@shu.edu.cn 







2 


Hengfei Ding, Changpin Li 


the process of mathematical modeling in the fractional realms, Caputo deriva¬ 
tives and Riemann-Liouville derivatives are mostly used. Generally speaking, 
the formers are often utilized to characterize history dependence, whilst the 
latter to describe long-range interactions. In contrast with the classical dif¬ 
fusion operator A, Riesz derivative operator, a special linear combination of 
the left Riemann-Liouville derivative operator and the right Riemann-Liouville 
derivative one, is applied to reflecting anomalous diffusion in space [l 7 [. The 

d a u(x) 


ath-order (1 < a < 2 ) Riesz derivative 
example, in [12; 


d\t 


in a; £ (a, b) is defined, for 


d a u(x) 

d\x\ a 


C a (rlD2,x + rlD%) u(x), 


( 1 ) 


where coefficient C a = — -— 4 -^—r, rlD% x and rlD x b are the left and right 

A COS l " 2 " O! I ’ 5 

Riemann-Liouville derivatives of order a defined by m 


1 


w(s)ds 


RlD°u(x) = < 


and 


r{2 — a) dx 2 J a (x — s) c 

d 2 u{x) _ 
dx 2 ’ ' 

1 d 2 f b u(s)ds 


—, 1 < a < 2, 


RLD X}b u{x) = < 


r(2 - a) dx 2 J x (s — x) a_1 ’ 
d 2 u(x) 

~d^’ a=2 ' 


1 < a < 2, 


The special case with a = —00 or b = +00 corresponds to the Liouville deriva¬ 
tive. For a well defined function on a bounded interval (a, b), we discuss them 
in [a, +00) or (—00,6] often by zero extension under suitable smooth condi¬ 
tions, i.e., let u{x) = 0 for all x > b or x < a. In this situation we have 
RlD* x u(x) = rlDZ^x) and RL D° b u(x) = r L D° +00 u(x). 

It is known that the Fourier transform of a given function u(x ) £ Li(IR) is 
given by, for example, in [Sj 


u(s) = T{u(x ); s} = 


p+oo 


e lsx u(a:)da;, x £ R, 


it follows that 


T 


rdjtW. ] 

\ dx n ’ J 


(is) n u(s), n £ IN, s £ R, 


( 2 ) 
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and 


T 



|s| a w(s), 1 < a < 2, seE. 


( 3 ) 


Note that —|s|“ = — (s 2 )2 for s £ R. So sometimes the Riesz derivative is 


,2 

also rewritten as a power of the operator —4^, i.e., 



Hence the Riesz derivative is often regarded as the symmetric fractional gen¬ 
eralization of the second derivative m- 


From (2) and (3), one easily sees that in the case a = 1, T i d Q\x\ I s \ ^ 



T | du ^ ; a|. Besides, Feller proposed another Riesz-type derivative (more 

general than Riesz derivative) with following form [2, 



with 



Letting the skewness parameter 0 = 0, one gets 



which is just the Riesz derivative (1). 

For most fractional differential equations, to obtain the analytical solutions 
are not easy even impossible, so many researchers have to solve fractional dif¬ 
ferential equations by using various kinds of numerical methods 



[28l[29lf30lf3Tll32 Il33ll34] . In particular, as for Riesz spatial fractional differential 


equations, the key issue is how to approximate the Riesz derivatives. From (1), 
one can see that a specific linear combination of the left and right Riemann- 
Liouvillc derivatives gives a Riesz derivative. So this question eventually come 
down to numerically approximate the Riemann-Liouville derivatives. Usually, 
we approximate the left Riemann-Liouville derivative by using the following 
Grlinwald-Letnikov formula 



due to the fact that Riemann-Liouville derivative and Griinwald-Letnikov one 
are equivalent under some smooth conditions m- But in specific applications, 
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we cannot solve a numerical problem with an infinite number of grid points, 
so one has to use the following formula 


hlDZMx) = — ™[ a eu{x - ih) + 0(h ), ( 4 ) 

£=0 

in which the Griinwald-Letnikov coefficients are given by 

_ (—iY ( a '\ — (—1Y _ + 1) _ p _ q i 

lj r(e+i)r(a-e+iy 

In fact, the generating function of the above coefficients vj^} is W\(z) = 

(!-«)“, i-e-, 

OO 

Wi(z) = (i - z) a = \z\ < 1. 

e=o 

Such coefficients can be recursively evaluated by 

1 + ol \ (a) 


(a) -1 (a) l 1 

®1,0 = 1, ^1,/ = I 1 - 


w 


= 0 , 1 ,... 


Unfortunately, it turns out to be unstable for the difference scheme for the 
time dependent equations by using (4) to approximate the Riemann-Liouville 
derivatives (or Riesz derivatives). In order to construct stable numerical schemes, 
one often needs to replace u(x — th) in (4) by u(x — (£ — p)h), where p £ 1R, 


rlD“ x u{x) = ™[“tu(x-(i-p)h) + 0{h), pY p ( 5 ) 


[^+ P ] 


r=o 


and 


, [^+p] 

RlD*, x u(x) = — w he u ( x ~ (£-p) h ) + 0{h 2 ), p = |, ( 6 ) 

t =o 

which is called as the shifted Griinwald-Letnikov formulas [T6j . 

At first sight one can find that the formula (6) has second-order accuracy. 
However, it needs some function values on nongrid points for the case a £ (0, 2) 
due to l — p IN. For the convenience of calculation and in order to avoid the 
nongrid point values by using the interpolation method, the optimal choose 
for p is: taking p = 0 for a £ (0,1] and taking p = 1 for a £ (1, 2 ). At this 
case, the shifted Griinwald-Letnikov formula (5) is used which gives lst-order 
accuracy. 

By combining the above shifted Griinwald-Letnikov formula, Tian et al. 
m developed two kinds of 2nd-order numerical schemes for the left Riemann- 
Liouville derivative as follows, 

RL D a,x U (x) = — 9 { ^}u{x -{e- 1 )h) + 0(h 2 ) 

f=0 
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and 


RR + 1 


RlDq X u(x) = Y g ( 2 ^u{x-{£-l)h) + 0(h 2 ), 


e=o 


where the coefficients </(') and g^ are given by 



(a) 

9i,o = 

a (a) 

2 ro i,o > 

(a) 

9i,l = 

a (a) 2 - a ( a ) 

2 ®i,< + 2 ,/-i» 

and 

(a) 

2 + a ( a ) 

(a) 

2 + a 

(a) (a) 2 + a ( a ) 

g 2 ,o - 

4 ro i, 0 i 

92,1 ~ 

4 " 

7 i,i i 9 2,1 ~ 4 w i,i 


> 1 , 


2-a 


,(«) 




i>2. 


On the other hand, the p-th order (p < 6) Lubich numerical differential 
formula 


RL 


d: 


, L ^ J 

ix u ( x ) = Y, w $ u ( x ~ ih )+ °( hP )> 


(7) 


e=o 


is derived by using the generating function below BE 


w P (z)=[Yi ( 1 -*) 1 


It should be pointed out that (7) holds for homogeneous initial conditions. 
The coefficients Y'l satisfy the following equation, 

/ P 1 \ a oo 

w p( z )= = X w p“M \ z \< 1 - 

\t=l / t=0 

The application of (7) to the spatial fractional differential equations with 
the Riemann-Liouvillc derivatives (or Riesz derivatives) is also unstable for 
a £ (1,2). To overcome this, we can propose the following shifted Lubich’s 
numerical differential formula, 

[^R 1 

RLDa, x u{x) = — Y w $ u ( x ~ (Z- l)h) + 0(h), p= 1 , 2 ,..., 6 . 

£=0 

But they have only lst-order accuracy by simple calculations. 

Because of the nonlocal properties of fractional operators, high-order nu¬ 
merical differential formulas lead to almost the same structure of the difference 
schemes as that produced by the lst-order scheme, but the former can greatly 
improve the computational accuracy. So it is more and more important and 
imperative to construct some effective and stable high-order numerical ap¬ 
proximate formulas. At present, the high-order numerical schemes are usually 
obtained by weighting the shifted and non-shifted Griinwald-Letnikov or Lu¬ 
bich difference operators In the present paper, our main goal is to 
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construct a class of much higher-order numerical differential formulas for Riesz 
derivatives by using another strategy. The key issue of the method is how to 
find the new class of the generating functions. The novelty of the paper is 
firstly to propose a 2nd-order formula for the Riemann-Liouville (or Riesz) 
derivatives based on its corresponding generating function, then developed 
the recurrence relations of the new generating functions. The main advantage 
of the method is the one can easily get unconditionally stable finite difference 
scheme. 

The paper is organized as follows. In Section 2, we derive a 2nd-order and 
several kinds of much higher-order numerical differential formulas for Riesz 
derivatives. In the meantime, the properties of coefficients, together with the 
convergence-order analysis of the 2nd-order formula are also studied. In Section 
3, the derived 2nd-order formula is applied to solve the Riesz spatial fractional 
advection diffusion equation. The solvability, stability and convergence anal¬ 
yses of the finite difference scheme are studied. Some numerical results are 
given in Section 4 in order to confirm the theoretical analyses. We conclude 
the paper with some remarks in the last section. 


2 New numerical differential formulas for Riesz derivatives 

In this section, we firstly develop a 2nd-order numerical differential formula 
for Riemann-Liouville derivatives and Riesz derivatives by using a new gen¬ 
erating function. Next, the properties of the 2nd-order coefficients have been 
discussed in details. Finally, the general forms of the much higher-order nu¬ 
merical differential formulas are also proposed. 

Theorem 1 Suppose u[x) £ C[“1+ 3 (]R,) and all the derivatives of u{x) up to 
order [a] + 4 belong to Li(R). Let 



Then if a = — oo, one has 


R L D-oo !X u(x) = L Bfu(x) + <D(h 2 ) 


(9) 


as h —> 0. 

Here [l = 0,1,...,) are the coefficients of the novel generating func- 




(10) 
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Proof Taking the Fourier transform on both sides of equation ( 8 ) yields 


where 


T{ l B%u{x)]s} = — 


{a) -i{t-l)hs 


u(s) 


fc =0 


= y—e ih °U 

h a 


( s ) 

(=0 

= (is) a (/>(ihs)u(s) 


(“) e -i ths 


e z ~ , _ 2a 2 — 6 a + 3 2 


<Kz) = = 1 - 


z 2 + 0(\zn 


So there exists a constant c\ > 0 satisfying 

\<j)(ihs) — 1 | < c\\s\ 2 h 2 . 


Furthermore, 

F{ l B 2 u(x)\ s} = (i s) a u(s) + (is) Q [</>(ihs) — 1 ]m(s) 

= T{ rlD^ ^x)-s} + <p{h, s), 

where <p(h,s) = (is)“[</>(ihs) — l]u(s). It follows that 

|£(M)| < ci|s| a+ 2 h 2 |u(s)|. 

Note that u(x) € C'[“1 + 3 (1R) and all the derivatives of u(x) up to order 
[a] + 4 belong to Li(lR). So there exists a positive constant c 2 such that 

|u(s)| < c 2 (l + |s|)- (W+4) . 

Taking the inverse Fourier transform of <p(h, s) yields 


\<fi(h,s) | = 


— J e' sx ip(h,s)ds 




< 


ClC 2 
27r 


(/ (l + |s|)“- [ “]- 2 ds^ h 2 = ch 2 , 


in which c = — jt—. --. Using again the inverse Fourier transform to 

7 r(|aJ + 1 — a) 

equation (11) gives 

RLD^ iX u(x)= L B%u(x)+0(h 2 ). 

This finishes the proof. 

By almost the same reasoning, one has the following theorem. 
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Theorem 2 Suppose u(x) £ C | [“]+”+ 1 (R) and all the derivatives of u(x) up 
to order [a] + n + 2 belong to L\ (ft). Then 


n— 1 


B?u(x) = rl D1^ x u{x) + Y, ( 7 £ “ rlD^M^)) h* + 0(h n ), n> 2. 


e =i 


e — 

Here f/ie coefficients 7 “ {£ = 1 , 2 ,...) satisfy equation —— W 2 (e~ z ) 

OO 

7 fz 1 , in which the coefficients of the first three terms are: 7 “ = 

j?=i 


2a 2 — 6 a + 3 
6 a 


3a 3 - 11a 2 + 12a-4 

12a2 


= 1 + 


0,7 


Ot 

2 


Next, we determine the coefficients of equation (10) by using the 
similar method presented in El- 



Comparing this equation with equation (10), one gets 
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With the help of equation (12) and automatic differentiation techniques [221 . 
one has the following recursive relations, 


,(«) 


3a- 2 
2a 


(q) _ 4a(l - a) j a) 
2 ’ x “ 3a-2 K2 ’ 0 ’ 


K S = ^ 3q 4 2 ) [ 4 (! - «)(“ - 1 


1 


(a) 


+(a-2)(2a-£ + 2)4J_ 2 


> 2 . 


(13) 


The above method is intuitive. Besides this, we can use another method to 
determine the coefficients k;") . Substituting 2 = e~ lx into (10), the coefficients 

K 2 I can be represented by the following integral form with the help of the 
inverse Fourier transform, 

1 _ 

/4°£ = 7T~ W 2 (-ix)e l(x dx, 

^ 7rl Jo 

where W 2 (-ix) = (^a=2 _ 2 (^i) e _ ix + £*=3 e -2i*)“ This type of integrals 
can be computed by the fast Fourier transform method Eli- 

Next, we study the properties of the coefficients (^ = 0,1,...). 


Theorem 3 The coefficients (£ = 0,1, ...) have the following properties 
for 1 < a < 2, 

"3a — 2' 


(i) 4“o = 


2a 


(a) _ 4a(l - a) (a) 

U, tv 2 i — ^ _ 2 ^ 2,0 ^ U ’ 


(ii) a 4 Q 2 = - (3a°- 2)^ 6a — ^ < 0 if a € (l,a*), while 


41 > 0 if a & [a*, 2), where a* = \ + ^S 21 + 48^ 

’ 8 24 

1.5333; 

(iii) k2°£ > 0 if £ > 3; 

(iv) 4“>~ J“Wr(° + i) ,-, „, 

’ 7T 

(v) 4°e —>■ 0 as £ — > 00 ; 


19 


y /621 + 48^87 


(vi) y = 0- 


£=o 


Proof (i) The direct computations give these results by formula (12). 
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(ii) With the help of the exact roots formula of cubic equation, one can 
easily get the conclusion. 

(iii) When l = 3,4, 5, we have the following results in view of (12), 

( a ) _ 4q(2 - a) (a - l) 2 /ii(a) ( a ) ( a ) a (a - l)(a - 2)fi 2 (a) („) 

' t2 ’ 3_ 3(3a — 2) 3 k 2,o>« 2,4- 6(3a — 2) 4 


and 


where 


j a ) _ 2a(2 - a) (a - l) 2 fj, 3 (a) (a) 


15(3a-2) 5 


v 2,0 > 


/ii (a) = 8a 2 — 7a + 2, 

/j 2 (a) = 64q 5 - 304q 4 + 507q 3 - 394q 2 + 148q - 24, 
fx 3 (a) = 64q 6 - 464q 5 + 1239q 4 - 1536q 3 + 984q 2 - 320q + 48. 

By simple computations, one has 

Hi(a) = 8q(q — 1) + a + 2 > 0, 

112 (a) = (a- l) 2 [q(8q - ll) 2 - 30q - 36] - 15(q - 1) - 3 < 0, 

113 (a) = (q - l) 2 (q - 2) [q(8q - 13) 2 - 82q - 20] + 53(q - l) 2 + 60(q - 1) + 15 > 0. 

So, K 2 ° 3 , K 24 an d k 25 are a d positive for 1 < q < 2. If £ > 6, we know that 
> 0 by the recurrence relation (13). It immediately follows that > 0 
for £ > 3. 

(iv) Using 


(- 1 )* 


sin (7rq) 


r(a - k + 1) 


r(k — a), 


the coefficients can be rewritten as 


( a ) sin (?rq) U(q + 1) (3a —2 
" 2 - < = n 


E 


q — 2 
3a- 2 


m— 0 

It is known that the ratio expansion of two gamma function 
N 


(a) r(£ — m — q) 
1 ’ rn r(£-m+ 1 )' 


r ( Z + Q ) _ a-b 

r(z + b) " 




Lfc=o 


holds as z —> 00 with | arg (z + a)| < 7r [26]. Here B^\a) are the generalized 
Bernoulli polynomials defined by lT8j 

U^T e “ = E^7 S i ff) («)’ 4 CT) («) = 1 ; \z\<2n, 
k / k =0 
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where (a) has the following explicit formula [23] 


4 a) (“) = E 


k\ (<7 + £ — 1\ l\ 


£—0 


xF 


W u 




f e (a+j) 


k-l 


£-k,£-a: 2£+l: 


a + j 


in which F[a, b;c;z] is the Gaussian hypergeometric function defined in [2] 

„ r , ab z a(a + 1)6(6 + 1) z 2 

= 1 + — yy + ^ ' -g*... 


So one has 


( a ) sin (7ra) F(a + 1) f 3a — 2 
K ‘ 2J - = n \~2^~ 


E 

m =0 


a- 2 
3a- 2 


,(«) 




Lfc=0 


Noting that 



as £ —> oo, 


one can get the coefficient k 


(a) 

2,1 


follows the power-law asymptotics, 


K 


(a) 

2,1 


sin (tto) T(a + 1) g _ a _ 1 
7r 


as £ —»■ oo. 


(v) From (iv), the asymptotics of k^} holds. Here we would rather use 

another approach to show it, where one can see that the Kifl is bounded by 

„(«) 
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One can see that 


(a) (a) 

W l,rn™lJ-r 


(a) ( 
W l,rn { 


X _ 1+0 


-m ) 1 


(a) 


» „(°0 


- 1±£L\ a ) 




> 1, m = 2,3,..., 


m +1 1 


Recalling that = 1) E°i = ~ a and > 0 for £ > 2 and 1 < a < 2, 
one gets 

(cc) (ck) ^ (ck) (a) 0 o 

< ro l,2 W i,/-2> m = 2,3, , 2 


It immediately follows that 


l 2/ 


< 


/3a- 2 

V 2a 


1 + 


2 — a 
3a-2 


(a) 

,/ - 


2 — a 


+ 


3a - 2 V 3a - 2 


2 — a 


i-1 


E 

m=2 

3a-2 


/ 2 — , 


(a) (a) 

nj, 0 tt7i 


2a 


\ 3a — 2 / 

) ot 

M(£,a)w[ a ^, £>2, 


where 

M(£,a) = 


1 + 


2 — a 
3a-2 


2 — a (2 — a 


a(3a — 2) / 2 — a 


3a-2 \3a — 2 


1 — 1 — a 


\3a-2j (£ — 1 — a)(£ — 2 — a) 


Because 


lim M(£,a) = l+ Q(2 a)(1 ° Q) >0, 1 < a < 2, 

^—>•+00 o(oQ: — 2 ) 

there exists a positive constant M(a), subject to M(£, a) < M(a) for 1 < a < 
2. In other words, we have 


» 

v 2,£ 


< M (a 




\ 2 a 


i,l ‘ 


So the 2nd-order coefficient r},fl is bounded by the lst-order coefficient • 

OO . . 

It is known that the positive series ^ w \t. is convergent [T3J- Therefore 

i=2 


the series K 2 *e i s also convergent. So the asymptotics of k'^ holds. 

3 =2 

(vi) By almost the same method used in na, the equality holds. 


» 
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Remark 1. For the right Liouville derivative, the approximation 
RL D« +00 u(x ) = R B^u(x) + 0(h 2 ), 

holds under the conditions Theorem ]! J Here, right difference operator R B 2 is 
defined by 

oo 

R B%u(x) = — {x + {£~ 1 )h). 


Remark 2. If u(x) is defined on [a, b\ satisfying the homogeneous condi¬ 
tions u(a ) = u(b) = 0, by suitable smooth extension one can get 

RlD% x u{x) = RLDfec x u{x) = L B 2 u{x ) + 0(h 2 ) = L A 2 u{x) + 0(h 2 ), 

(14) 

and 


RL 


D x , b u ( x ) = RL D^ +00 u(x) = R B«u(x) + 0(h 2 ) = R A a 2 u{x) + 0]h 2 ). 


Here the operators L A 2 and R A 2 are defined as follows, 


(15) 


[*7^1 + 1 


L A 2 u{x) = — K i a e u l ) h ) 


t =o 


and 


r A 2 u{x) = i ^ k^Ju {x + {t- 1 )h). 


e=o 


Hence, combining equations (1), (14) and (15), one can obtain a new kind 
of 2nd-order difference scheme for Riesz derivatives (1), 

= C <* ( + R A 2 u{x)) + 0(h 2 ). ( 16 ) 

Finally, we give the more general high-order numerical algorithms below. 

Theorem 4 Let u(x) £ C'[“1 +P+1 (R) and all the derivatives of u{x) up to 
order [a] + p + 2 belong to Li(lR). Set 

1 OO 

l B*u{x) = — K< p a e u 1 )h ), 

fcO 

and 

1 OO 

R BpU(x) = --J2 Kpju {x + {t- 1 )h). 

fco 


Then 


RlD- oo, x u (x) = L Bpu(x) + 0(h p ), p> 3, 
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and 


rl D^ +00 u{x) = R BpU(x) + 0{hP), P > 3. 


Here the generating functions with coefficients (l = 0,1,...) for p > 3 


are 


W P {z) =(!-*) + 


a —2 
2 a 


a-~E+E 


v(«) 

A /c-l ,k—l 


k =3 


a 


(i-^ 


i.e., 

oo 

Wp(^) = E N < !> P > 3. 

£=o 

in which the parameters AjE (fc = 3,4,...) can 6e determined by the 
following equation 


W k , s (e 


*) — = 1 - 

' z a 


OO 


\(“V 

/ , a M " 
t=k 


k = 2,3,... 


Proof The proof of this theorem is almost the same as that of Theorem [TJ so 
we omit it here. 

Similarly, define the following p-th order difference operators 


[^rd+i 


L Apu(x) = — E K p} u i x -(*- l ) h ) 


and 


£=o 


[Vl+i 


R ApU(x) = E E K { p}u(x + {i-\ )h) 


i -o 


then the p-th order numerical differential algorithm for Riesz derivatives in 
(a, b) is given by 


+ *^«(*)) + <X P ), P > 3- 


The cases p = 3 and p = 4 are listed in Appendix A for reference. 
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3 Application of the 2nd-order scheme 


In this section, we apply the derived 2nd-order scheme to the Riesz space 
fractional partial differential equation. 

We study one-dimensional Riesz spatial fractional advection diffusion equa¬ 
tion in the following form, 


du(x,t) du(x,t) d a u(x ) 

dt + dx “ d\x\ a 


f{x,t), a < x < b, 0 <t <T, (17) 


with the initial condition 


u(x, 0) = u°(x ), a < x < b, 
and the Dirichlet boundary conditions 

u(a, t ) = u(b, t) = 0, 0 < t < T, 


where K > 0 and K a > 0 are the advection and diffusion coefficients, respec¬ 
tively. f(x,t) and u°(x) are suitably smooth functions. 

Let Xj = jh ( j = 0,1, ■ ■ • , M) and tk = hr (k = 0,1, • • ■ , N), where h = 
and t = jj are the uniform spatial and temporal meshsizes, respectively. 
And M, N are two positive integers. Denote it* = u{xj,tk), 0 < k < N, 0 < 
j < M, then the computational domain [0,T] x [a, b] is discretized by Q T h = 
fl T x flh, where 1? T = {tk | 0 < k < N} and fih = {xj\ 0 < j < M}. Given any 
grid function [uj\ 0 < j < M, 0 < k < N} on Q t h, denote 


k--k 


= 2 + U J 


k -1 




k -1 


), 


u = l ("? - . *«<•? = \ (4»‘ + . + «,»*_.) . Si „J = i (4»J + . - s,u)_ h 


For convenience, let Vh = {u| u = {itj| 0 < j < M} is a grid functions on 
fih and uq = um = 0}. Then for any grid function u, v £ Vh, we can define 
the following inner products 


M—l M 

(u,v) = h UjVj, (S x u,S x v) = (tixUj- i) (s x Vj_^) , 


3 =1 


and the corresponding norms 


3=1 


IMI = t/{u,u), ||4m|| = y/(S x U,S x u). 

Next, considering equation (17) at the gird points (xj , t k+ \), one has 
du{xj,t k+ 1) du(xj,t k+ 1) a Q u(a; J ,t fc+ i) 


dt 


dx 


= K n 


d\x\ 


+ f(Xj,t k+ i). 
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Substituting (16) into the above equation leads to 
du(xj,t k+ 1) 


at 


+ KS x u(xj, t k+ 1) — K a S x ii(xj, ) + f (xj + 0(h ), 


where the operator 5% is defined by 6 %u(x,t) = C a ( L A% + u{x,t)- 

Using Taylor expansion yields 


du(xj,t k + 1 ) 
at 


= S t u(x j ,t k +i) + 0(T 2 ). 


(18) 


A combination of the above two equations gives, 

5tu{xj, tfc-|-i) + K5xu(xj, tfe_|-i) = K a 5 x u(xj, t k+ 1 ) + f(xj,tk+x)^~Rj> 

where there exists a positive constant C3 such that 

\Rj \ < c 3 (r 2 + h 2 ), 0 < k < N — 1, l < j < M — 1. 

Omitting the small terms in (18), and replacing the grid function 

fc+i 

u(xj,tk+ i) with its numerical approximation U 3 2 , we obtain the following 

finite difference scheme for equation (17), 


S t U- +i + K6 S U- + ^ = K a 5«U k 3 +h + f- +i 
k = 0,1,..., AT — l,j = 1,2, ... ,M — 1, 


(19) 


U° = u°(xj), j = 0,1,...,M, 

U^ = U^ = 0, A; = 0,1,...,TV. 

We now prove the solvability, stability, and convergence of the difference 
scheme (19). Firstly, let us list some preliminary results. 

Definition 1 |Tj Let n x n Toeplitz matrix T„ be in the form: 

/ to f_i ••• t2-n tl-n\ 
tl to t -1 ■ ■ ■ t2-r 


T.„ = 


: ti to '■ : 

t n -2 ••• '■ '■ t-i 

tn-l tM -2 • • • tl to 


i.e., t,; ; = ti-j. Assume that the diagonals are the Fourier coefficients of func¬ 
tion /, i.e., 

tk = 2tt / ^( a; ) e_lfcXdx ’ 

then function /(a;) is called the generating function of T„. 
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Lemma 1 (Grenander-Szego Theorem fffj) For the above Toeplitz matrix T n , 
let f(x) be a 2 tt- periodic continuous real-valued function defined on [—7r, 7r]. 
Denote A m i„(T„) and A max (T„) as the smallest and largest eigenvalues of 
T„, respectively. Then one has 

fmin A m i ra (T„) ^ XmaxifFn) — fmax, 

where f m in, fmax are the minimum and maximum values of f{x) on [—7r,7r]. 
Moreover, if fmin < fmax, then all eigenvalues of T„ satisfy 

fmin ^ A(T 

n) ^ fmax, 

for all n > 0. And furthermore if fmax < 0? then T n is negative semi-definite. 


Theorem 5 Denote 


G 


a 


(a) 

k 2,1 

(a) 

^2,0 

0 


0 

(a) 

k 2,2 

(a) 

k 2,1 

(a) 

k 2,0 

0 


(a) (a) 

'2,M—2 K 2,M-3 


(a) 

*2,1 

(a) 

K 2,0 

.(a) , («) 

2.M—1 K 2,M-2 


(a) 

k 2,2 

(a) 

K 2,l 


Then matrix G = (G Q + G„) is negative semi-definite. 

Proof According to Definition [I] we know that the generating functions of the 
matrices G a and G^ are 


f Ga (x) = J2 K{ 3 ei{e ~ 1)x ’ /gz(*) = £' 


£=0 


"2 ’ 


£=0 


respectively. Accordingly, the generating function of matrix G is f(a,x) = 
/g q ( a )+/ g t (a;), which is a periodic continuous real-valued function on [—7r, 7r]. 
Application of equation (10) leads to 


f(a,x) = e~ 


3 a 2 2 (a l) ^ _|_ ®_ - e 2ix 


2 a 


a 


2 a 


+e 


3® 2 2(q 1) £~ix _j_ ® _ -c~ 2ix 


2 a 


a 


2 a 



— — -e ix ] 
3a- 2 J 


+e ix (1 



-— -e^A 

3a- 2 J 


Since /(a, x) is an even function, we only need consider its principal value 
on [0,7r]. Using the following formulas 
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and 


(l_ e ±ix)“ = ^2sin|) a e ±ia (^) 


(a — bi) a = (a 2 + 6 2 ) 2 e la6 , 9 = — arctan —, 


then one gets 

,, . / . x\ a ( 3a — 2 

f{a,x) = (2sm-J 

where 


2a 


a — 2 

1-COS X 

3a- 2 


a — 2 
3a- 2 


■ sin x 


Q(a, x) = 2 cos ( a ( 9 + 


and 


Let 


Then 


(ot — 2) sin x 

11 = -‘‘ rCt ‘‘ n (3a~2) — („ — 2)cos* ' “ £ (1 ’ 2 >' X 6 [M ' 


Z(a, x) = a ( 9 + ——— I — x, at (1,2), x £ [0,7r]. 


d0 a 1 _ 2(1 — a)(2 — a)(3a — 2) sin 2 (|) 

X (a,x) + 2 4(a — l) 2 + (2 — a)(3a — 2) cos 2 (§) — ’ 

that is to say that Z x (a. x) is an monotonically nonincreasing function with 
respect to x, so 


Z max (a,x) = Z(a, 0) = - -a, 


and 


Qmax(d, 3?) — 2 COS (^max(^) ^) ) — 2 COS ^ ^ a^ < 0, a (1,2). 


Hence, we know that f(a,x) < 0 and matrix G = (G a + Gf) is negative 
semi-definite for a £ (1, 2) by Lemma[TJ 

Theorem 6 For any v £ Vh, the following inequality holds 

(6%v,v) <0 for a£ (1,2). 


Proof One can easily check that 

{6y,v)=C a (( L A%+ R Af) v,v) = C a hv T (G a + Gl)v = C a hv T Gv, 
which implies that (6%v,v) < 0 by Theorem [5] 

Theorem 7 Finite difference scheme (19) is uniquely solvable for 1 < a < 2. 


Q(a,x), 
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Proof Here we use induction method to show it. From (19), it is obviously 
that the result holds for k = 0. 

Now suppose that U k has been determined by equation (19) for 1 < k < 
N- 1, i.e., 

6 t U- +i + K5 s U k 3 +h = K a 5<fU k+h + 

which can be rewritten as 

2U k+1 + rK5xUj +1 - TK a 5fUj +1 = 2U k - TK5 s U k + TK a 5*U k + 2 rf k+ K 

Considering the homogeneous form of the above equation and taking the 
inner product with U k+l yield 

2(U k+1 ,U k+1 ) + rK(6xU k+1 ,U k+1 ) - rK a {6^U k+1 ,U k+1 ) = 0. 


Because 


(6xU k+1 ,U k+1 ) =hJ2 (W* +1 ) U k+1 = - Y, - U j-i) U j +1 = 0- 

i=i 1 l=i 

and 

(5fU k+1 ,U k+1 ) < 0, 

we have 

(U k+1 ,U k+1 ) < 0. 

So, ||t/ A;+1 || =0 and U k+1 can be solved uniquely. 

Theorem 8 Finite difference scheme (19) is unconditionally stable with re¬ 
spect to the initial values for 1 < a < 2. 

Proof Suppose that v k is the solution of the following difference equation, 


Stv J +l + KSxV k+i = K a sy k+i + f 3 +i , 
k = 0,l,...,N—l,j = l,2,...,M—l, 


( 20 ) 


Vj =u°{xj) + p j7 j = 0,1,..., M, 
v o= v m=°i k = 0,l,...,N. 

Let f k = U k — v k , then from equations (19) and (20) one has 

5ff k+h + KSxC" 2 = 

k = 0,1,... ,N — l,j = 1,2, ..., M — 1, 

$ = P„ j = 0,l,...,M, 

Co=Cm = 0, k = 0,l,...,N. 


(21) 


Denote 


£ = (£i> ■ • ■ i Cm-i), P = (pi, ■ ■ ■, Pm- i). 
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Taking the inner product of (21) with £ fe+ 2 yields 

(<54 fc+l ,€ fe+ ^) +if(4£ fc+l ,€ fc+l ) =K a (s°e k+ t,e h+ t). 

Note that 

( 5t € k+i ,£ k+i ) = (€ fe+1 -£ fc ,£ fe+1 + £*) = ^ (ll£ fc+1 ll 2 - ll£ fc ll 2 ) 

M—1 M—1 

(<k£ fc+ *,£ fc+ *) =/i E (5^ +i )^ +4 

and 


1=i 


1=1 

||^ + 1 || 2 -||^|| 2 <o, 

||€ fc+1 ||<||€°|| = ||p|l> 

that is to say that finite difference scheme (19) is unconditionally stable with 
respect to the initial values. All this finishes the proof. 

Theorem 9 Finite difference scheme (19) is convergent with order 0 (t 2 + 
h 2 ). 


then we have 


i.e., 


Proof Suppose that u(xj,tk ) be the exact solution of equation (17) and U k 
be the solution of difference equation (19). Let e k = u(xj,tk) — U k , then from 
equations (17) and (19), one gets 


5 t e) +h + K5= K a 5fe) + ^ + R k , 

k = 0,1,...,1V — l,j = 1,2,..., M — 1, 

e? = 0, j = 0,1,..., M, 
e k 0 =e k M = 0, k = 0,1,... ,N. 


Set 


e = (ei,..., £m- i); R = (Ri, ■ ■ ■ j Rm- i)- 


Taking the inner product of (22) with e k+ l leads to 


( 22 ) 


((5 t e fc+ 5, £ fc+ ^ + K (d s e k+ i,£ k+ ^ = K a ^“e fe+ 5, £ /c +2^ + (R k , £ k +h^ . 

(23) 

Since 


(R k ,e k+ i) < ||i? fc || 

e k +h 

0 $ 

II 

£ k + £ k +1 

V ’ ) w w 



2 


(24) 
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we have the following estimate in view of (23) and (24), 


id 


£ fc+1 || - 11 £T fc 


04 


lie’ll ll £ fc + £ fc + 1 l 


i.e., 


|e fe+1 1| < ||e fc || + r ||i? fc || < ||e°|| + r ^ ||l? n || . 


n —0 


Notice that 


M—1 


|l ? n || 2 = hJ 2 (Rj) < (M - 1 )hcl(T 2 + h 2 ) 2 <{b- a)c 2 (r 2 + h 2 ) 2 . 
3 = i 


Then 


r ^ ||i?"|| < kry/b — clcz{t 2 + h 2 ) < c^Ty/b — a(r 2 + h 2 ), 


n —0 


which gives 

||e fc || < Ci(r 2 + h 2 ), \ < k < N, 
where C4 = c^Ty/b — a. This ends the proof. 


4 Numerical examples 


In this section, we present one numerical example for checking the conver¬ 
gence order of the numerical formula (16). Next, another numerical example is 
given to test the convergence order and numerical stability for finite difference 
scheme (19). 


Example 1 Consider function u(x) = x 2 (l — x) 2 , x £ [0,1]. The Riesz deriva¬ 
tive of u(x) at x = 0.5 is 


d a u(x) ( a 2 — 6a + 8)2“ 1 

d\x\ a ^ =0 ' 5 = r{ 5- a) 



Table 1 lists the absolute errors and numerical convergence orders at x = 
0.5 with different a and step size h. From the results presented in Table 1, one 
can see that the convergence orders are in line with the theoretical analysis. 


Example 2 We consider the following Riesz spatial fractional advection-diffusion 
equation in the following form, 


1 du(x,t) n du(x,t ) 2 <9“u(x) 

dt + dx =a d\x\ a + 


f(x,t ), 


0 < x < 1, 0 < i < 1, 
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Table 1 The absolute errors and convergence orders of Example 4.1 by numerical scheme 
(16). 


Q 

h 

the absolute errors 

the convergence orders 

i.i 

l 

20 

2.492284e-003 

— 


1 

40 

6.462793e-004 

1.9472 


1 

80 

1.643881e-004 

1.9751 


1 

160 

4.144518e-005 

1.9878 


1 

320 

1.040456e-005 

1.9940 

1.3 

1 

20 

3.563949e-003 

— 


1 

40 

9.146722e-004 

1.9621 


1 

80 

2.315235e-004 

1.9821 


1 

160 

5.823146e-005 

1.9913 


1 

320 

1.460130e-005 

1.9957 

1.5 

1 

20 

4.555022e-003 

— 


1 

40 

1.157683e-003 

1.9762 


1 

80 

2.916709e-004 

1.9888 


1 

160 

7.319217e-005 

1.9946 


1 

320 

1.833193e-005 

1.9973 

1.7 

1 

20 

5.266851e-003 

— 


1 

40 

1.326934e-003 

1.9888 


1 

80 

3.329312e-004 

1.9948 


1 

160 

8.337785e-005 

1.9975 


1 

320 

2.086231e-005 

1.9988 

1.9 

1 

20 

5.352412e-003 

— 


1 

40 

1.339793e-003 

1.9982 


1 

80 

3.351414e-004 

1.9992 


1 

160 

8.380846e-005 

1.9996 


1 

320 

2.095494e-005 

1.9998 


with the source term 
f(x,t)=a 2 ^ 12 

+ 


r(5 - a) 
2160 


[x 4 ~ a + (1 - x) 4 ~ a ] - 


240 


r(6- a) 


[x 5 ~ a + (1 - z) 5 “ a ] 


+ 


m-a) 

20160 


[x 6 ~ a + (1 - x) 6 ~ a ] - 


10080 
r(8 - a) 
2 


[x 7 ~ a + (1 - x) 7 ~ a ] 


lx 8 “ + (1 — a;) 8 a l 1 cos ^ | — 2 atx 4 (l — x) 4 sin (at 2 ) 
J cos (fa) 


r(9 - a) 

+2 cos(at 2 )(8x 7 — 28x 6 + 36x 5 — 20x 4 + 4* 3 ). 


Its exact solution u(x, t) = cos(at 2 )x 4 (l—x) 4 satisfies the corresponding initial 
and boundary values conditions. 
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In order to check the convergence order in temporal direction, we apply 
finite difference scheme (19) on a fixed sufficiently small spatial stepsize h and 
variable temporal stepsizes r. Similarly, in order to check the convergence or¬ 
der in spatial direction, we use a fixed sufficiently small temporal stepsize r 
and variable spatial stepsize h. The absolute errors and numerical convergence 
orders are presented in Tables 2 and 3, respectively. From these numerical re¬ 
sults, it is seen that numerical scheme (19) has 2nd-order convergence order 
for both temporal and spatial directions, which is in agreement with the de¬ 
rived theoretical results. Furthermore, in Figs. |T] and [5] we present the errors 
for different r, h and a, which also show that finite difference scheme (19) is 
effective. 


Table 2 The absolute errors and temporal convergence orders of Example 4.2 by numerical 
scheme (19) with h = x . 


a 

T 

the absolute errors 

the spatial convergence orders 

1.2 

1 

5 

8.853323e-005 

— 


1 

10 

2.139870e-005 

2.05 


1 

20 

5.374545e-006 

1.99 


1 

40 

1.350102e-006 

1.99 


1 

80 

3.461897e-007 

1.96 

1.4 

1 

5 

9.968682e-005 

— 


1 

10 

2.551591e-005 

1.97 


1 

20 

6.323861e-006 

2.01 


1 

40 

1.591389e-006 

1.99 


1 

80 

4.064288e-007 

1.97 

1.6 

1 

5 

1.238098e-004 

— 


1 

10 

2.974978e-005 

2.06 


1 

20 

7.370363e-006 

2.01 


1 

40 

1.846014e-006 

2.00 


1 

80 

4.695380e-007 

1.98 

1.8 

1 

5 

1.435874e-004 

— 


1 

10 

3.361038e-005 

2.09 


1 

20 

8.398948e-006 

2.00 


1 

40 

2.099550e-006 

2.00 


1 

80 

5.309125e-007 

1.98 
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Table 3 The absolute errors and spatial convergence orders of Example 4.2 by numerical 
scheme (19) with r = . 


a 

h 

the absolute errors 

the spatial convergence orders 

1.2 

l 

10 

2.101375e-004 

— 


l 

20 

5.260133e-005 

2.00 


1 

40 

1.334869e-005 

1.98 


1 

80 

3.373047e-006 

1.98 


1 

160 

8.484737e-007 

1.99 

1.4 

1 

10 

2.015275e-004 

— 


1 

20 

5.155201e-005 

1.97 


1 

40 

1.312357e-005 

1.97 


1 

80 

3.315909e-006 

1.97 


1 

160 

8.337871e-007 

1.99 

1.6 

1 

10 

1.831026e-004 

— 


1 

20 

4.672567e-005 

1.97 


1 

40 

1.183333e-005 

1.98 


1 

80 

2.979178e-006 

1.99 


1 

160 

7.475850e-007 

1.99 

1.8 

1 

10 

1.485772e-004 

— 


1 

20 

3.783943e-005 

1.97 


1 

40 

9.534022e-006 

1.99 


1 

80 

2.391298e-006 

2.00 


1 

160 

5.991744e-007 

2.00 


5 Conclusion 


In this paper, a class of new numerical schemes are proposed for Riemann- 
Liouville derivatives (and Riesz derivatives). Application of the 2nd-order 
scheme in spatial direction and the Crank-Nicolson technique in temporal 
direction, a finite difference scheme is developed to solve the Riesz space frac¬ 
tional advection diffusion equation. We prove that the difference scheme is un¬ 
conditionally stable and convergent by using the energy method. Finally, two 
numerical examples have been given to show the effectiveness of the numerical 
schemes. Following this idea, the method and technique can be extended in 
a straightforward way to construct much higher-order numerical algorithms 
for the tempered and substantial fractional derivatives, and the corresponding 
fractional differential equations. 
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—. 

i 


i 


Fig. 1 The error surface between the exact solution and numerical solution with r = -L 


Appendix A 

Now we present cases p = 3 and p = 4 in details. 

(i) P = 3 

The generating function with coefficients {t = 0,1,...,) reads as, 



W 3 (z) = (a 3 1 + a 32 z + a 33 z 2 + a 34 z 3 )“ = ^ k$ z e , 

1=0 


where 


031 = 


11a 2 — 12a + 3 
6a 2 

3a 2 — 8a + 3 
2a 2 ’ 


■ i ^32 = 


—6a 2 + 10a — 3 
2^ ’ 
-2a 2 + 6a - 3 
6a 2 


033 = 


034 = 
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Fig. 2 

and h = 


The error surface between the exact solution and numerical solution with r = 
i 

1000 • 


1 

100 


This generating function can be also rewritten as 
W 3 (z) = (a 3 1 + a 32 z + a 33 z 2 + a 34 z 3 ) a = a 31 (1 - z) a (l + 6312: + b 32 z 2 ) 


2\ a 


ii 


■ 031 (i-zfE (?) ( & 3i^) fi EE 

4=0 v iy 


=0 


332 

^31 


^2 


^31 


E 

1=0 


( 1 \ 4 “ t “^ 2 (Q [) \| 7 ^ 1 — 2^2 ^2 

V" (-1) (^1-^2)!o 3 i O 32 («) (a) 

2^ f n \(P, — o/L'ii w u-4 w i,4-4 


in which 


So, 


&31 = 


IT 1 ] 

E 

0 £2—0 

— 7a 2 + 18a — 6 


11a 2 - 12a+ 3 ’ 


[2^1] 


b 32 = 


2a 2 — 6a + 3 
11a 2 - 12a+ 3' 


44 = °3i E E £ = °> 1 .• • •• 

4=0 4=0 
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where 


-faf M (- 1)^(4 -4)! u,-v a&& 
P(a ’ £l ’ £2) “ 41(4-24)1 631 &32 ‘ 


In addition, we can get the following recursion relation by using the ex¬ 
pressions of Kg Q £ and automatic differentiation techniques, 

’ I la 2 - 12a+ 3' 


,(«) 


6 a 2 


t a ) 3a(6a 2 — 10a + 3) ( a \ 

K3 ’ 1 = 11a 2 -12a+ 3 K3 ’ 0 ’ 

(a) 3a(108a 5 - 402a 4 + 520a 3 - 312a 2 + 87a - 9) (a) 

K3 > 2 ~~ 2(lla 2 — 12a + 3) 2 Ks ’ 0 ’ 


(a) _ 

~ a 31 e L 


^32 (ct — £ 1)^3 d - 033(2a — £ -\- 2)^3 A—2 


(a) 


+034 (3a — £ + 3)K3^_3 


, l > 3. 


(ii) p = 4 

The generating function with coefficients {£ = 0,1,...,) reads as fol¬ 


lows, 


W±{z) = (041 + 0422 + O432 2 + 0442 s + 04 5 2 4 )“ = ^ K 4^2 £ , 


t =0 


in which 


25a 3 - 35a 2 + 15a - 2 -24a 3 + 52a 2 - 27a + 4 

«41 = -72 o-, O42 = - 


12 a 3 

6 a 3 - 19a 2 + 12a- 2 
043 = - ^ -, a 44 = 

3a 3 - 11a 2 + 9a-2 
“ 45 = - Utf -■ 

Similarly, one can get 

z [i c A [h E A 


6 a 3 

—8a 3 + 28a 2 — 21a + 4 
6 a 3 


l\ —0 ^ 2—0 £3=max{0,2£2-^i} 


,(«) 


where 


P(a,e i,4,4) 


(~l) <1+<a (4-4)! 

4! (4 -24)! (4 + 4 -24)! 


h tl+l3-2t, 2h t. 2 -2t 3h f. 3 
°41 °42 °43’ 


and 


641 

&43 


—23a 3 + 69a 2 — 39a + 6 
25a 3 - 35a 2 + 15a - 2 ’ 
—3a 3 T 11a 2 — 9a T 2 
25a 3 - 35a 2 + 15a- 2' 


13a 3 - 45a 2 + 33a - 6 
25a 3 - 35a 2 + 15a-2’ 
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The recursion formula is given as, 


' (a) _ /25a 3 - 35a 2 + 15a-2 


l 


12a 3 


( Q ) 2a(24a 3 — 52a 2 + 27a — 4) ( Q ) 


k 4,i = - 


25a 3 - 35a 2 + 15a - 2 4 ’° ’ 

(a) _ 2a(576a 7 - 2622a 6 + 4441a 5 - 3835a 4 + 1844a 3 - 497a 2 + 70a - 4) (o 


M.2 


(a) 

o = - 


(25a 3 - 35a 2 + 15a - 2 ) 2 


M,0 J 


2 a 


(27648a 11 - 19785a 10 + 591000a 9 - 995240a 8 


4 ’ 3 3(25a 3 - 35a 2 + 15a - 2 ) 3 

+1067901a 7 - 775354a 6 + 390051a 5 - 135738a 4 + 31923a 3 - 4820a 2 
+420a - 16) f4“o, 


( a ) _ 

~ a 4 !l L 


a 4 2 (a — l + + 043(20 — i + 2)k k ^J_ 2 + 044(30 — t + 3)k^_ 3 

> 4. 


» 


„(«) 


+045 (4a — £ + 4)k 


(a) 

4, e -4 


The cases for p > 5 can be similarly derived howbeit very complicated. We 
omit them here. 
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